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ABSTRACT 

We explore magnetohydrodynamic (MHD) solutions for envelope expansions with core 
collapse (EECC) with isothermal MHD shocks in a quasi-spherical symmetry and out- 
line potential astrophysical applications of such magnetized shock flows. By including a 
random magnetic field in a gas medium, we further extend the recent isothermal shock 
results of Bian & Lou who have unified earlier similarity isothermal shock solutions of 
Tsai & Hsu, of Shu et al. and of Shen & Lou in a more general framework. MHD shock 
^^ ' solutions are classified into three classes according to the downstream characteristics 

f*^ , near the core. Class I solutions are those characterized by free-fall collapses towards 

the core downstream of an MHD shock, while Class II solutions are those character- 
ized by Larson-Penston (LP) type near the core downstream of an MHD shock. Class 
HI solutions are novel, sharing both features of Class I and II solutions with the pres- 
Q ■ ence of a sufficiently strong magnetic field as a prerequisite. Various MHD processes 

H I may occur within the regime of these isothermal MHD shock similarity solutions, such 

C/3 i as sub-magnetosonic oscillations, free-fall core collapses, radial contractions and ex- 

J^ ' pansions. Both possibilities of perpendicular and oblique MHD shocks are analyzed. 

Under the current approximation of MHD EECC solutions, only perpendicular shocks 
, , are systematically calculated. These similarity MHD shocks propagate at either sub- 

^S ' magnetosonic or super-magnetosonic constant speeds. We can construct Class I, II 

and HI MHD shocks matching with an isothermal magnetostatic outer envelope or an 
MHD breeze. We can also construct families of twin MHD shock solutions as well as 
an 'isothermal MHD shock' separating two magnetofluid regions of two different yet 
constant temperatures. The versatile behaviours of such MHD shock solutions may be 
utilized to model a wide range of astrophysical problems, including star formation in 
magnetized molecular clouds, 'MHD champagne flows' in HII regions around luminous 
massive OB stars, MHD link between the asymptotic giant branch phase to the proto- 
planetary nebula phase with a hot central magnetized white dwarf, relativistic MHD 
pulsar winds in supernova remnants, radio afterglows of soft gamma-ray repeaters and 
so forth. 

Key words: magnetohydrodynamics (MHD) - radiation mechanisms: general - shock 
waves - stars: AGB and post-AGB - stars: formation - stars: winds, outflows 



1 INTRODUCTION gravitational inflows and outflows during a certain phase 

of evolution. These problems with different approximations 



S 



Many astrophysical processes, including star formation 
stellar collapse, supernova explosions, gamma-ray burst; 
(GRBs) and galaxy cluster evolution etc., may involve self- 



(e.g., isothermal or polytropic equations of state, spherical 
stellar collapse, supernova explosions, gamma-ray bursts i-i-i ^ ■ i-^- i-^- 

^ . ° . -^ or cyhndncal symmetries, radiative or non-radiative regimes 



etc.) have been explored extensively in various contexts (Bo- 
denheimer & Sweigart 1968; Shu 1977; Goldreich & Weber 
1980; Suto & Silk 1988; Foster & Chevalier 1993; Roily & 
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Shadmehri 2005). When a gas medium has evolved suffi- 
ciently far away from the influence of initial and bound- 
ary conditions, it may gradually evolve into self-similar be- 
haviours (e.g., Sedov 1959; Landau & Lifshitz 1959; Baren- 
blatt & Zel'dovich 1972) . This evolving trend as indicated by 
numerical calculations to adjust to a roughly self-similar pro- 
file led Larson ( 1969a, b) and Pension ( 1969a, b) to develop 
a spherically symmetric isothermal self-gravitating similar- 
ity solution which has a flat inner density profile with a 
radial speed proportional to radius and a supersonic outer 
envelope with an expansion speed 3.3 times the isothermal 
sound speed. Shu (1977) constructed an alternative self- 
similar expansion- wave collapse solution (EWCS), in which 
the initial condition is an isothermal sphere with a gas mass 



density p oc 



throughout in a hydrostatic equilibrium. 



Perturbations or some central energy loss initiate the col- 
lapse and a rarefacetion wave front travels outward through 
the gas medium inside of which gas materials rapidly at- 
tain free-fall velocity towards to the centre. This EWCS sce- 
nario of inside -out collapse for protos t ar form ation has been 
advocated by IShu. Adams fc Liza no' (1987) and compared 
with observations JA dams^ Lada & Shu 1987; Zhou et a l. 
ll993l:IChoi et al.ll995l:ISaito et al . 1999; H arvev et alj200ll) . 
Immediately following the analysis of Shu (1977), Hunter 
(1977) constructed complete isothermal self-similar solu- 
tions by going from negative times to positive times and pro- 
posed a new matching method to find solutions in density- 
speed phase diagrams. He found an infinite number of dis- 
crete solutions in the 'complete solution space'. Whitworth 
& Summers (1985) have further expanded these solutions 
into two-parameter continua by allowing weak discontinu- 
ities across the sonic critical line. Expansion wave solutions 
for more general adiabatic and self-gravitating gas fiows were 
discussed by Cheng (1978). Suto & Silk (1988) studied gener- 
alized poly tropic similarity solutions. Yahil (1983) obtained 
polytropic fiow solutions of the Larson-Penston (LP) type 
with the polytropic index 7 in the range of 6/5 < 7 < 4/3 
for applications to stellar collapse and noted that various 
solution features appear to be i n accordance with results 
of numerical simulations. iFoster fc Chevalier (1993) studied 
the gravitational collapse of an isothermal cloud by hydro- 
dynamic simulations. They recovered the LP solution in the 
central r egion where a core forms and the self-similar so- 
lution of IShul (|l973) when the ratio of initial outer cloud 
radius to core radius is > 20. 

Shock processes can naturally occur in diverse astro- 
physical settings, for example, supernova explosions, pho- 
toionized gas, stellar winds, collisions between high velocity 
clumps of interstellar gas, collisions of two or several galax- 
ies etc. Shocks have been under extensive investigations in 
contexts of supernova remnants (SNRs), nebulae associated 
with T Tauri stars, planetary nebulae, HII regions, parti- 
cle accelerations and GRBs etc. (e.g., McKee & HoUenbach 
1980; Draine & McKee 1993; Meszaros 2002). McKee & Os- 
triker (1977) explored the interaction of a SNR with an in- 
homogeneous ambient medium. Self-similar hydrodynamic 
shocks were studied for interaction zones of colliding stellar 
winds (Chevalier 1982; Chevalier & Imamura 1983). Tsai & 
Hsu (1995) constructed self-similar shocks describing a sit- 
uation in which a central thermal or kinetic energy release 
initiates an outgoing shock during a protostellar collapse to 
form a low-mass star. Tsai & Hsu (1995) also constructed a 



shock expansion solution with a finite core density matched 
with a static SIS envelope; this solution has been generalized 
into a family of 'champagne flow' shock solutions by Shu et 
al. (2002) to model expansions of HII regions surrounding 
massive OB stars after a rapid passage of an initial ioniza- 
tion front in a neutral hydrogen cloud. Shen & Lou (2004) 
studied EECC shock solutions by allowing outflows or in- 
flows far away from the central region. In terms of modelling 
protostellar systems, this flexibility can accommodate a va- 
riety of physical possibilities. Bian & Lou (2005) explored 
the parameter space more systematically to construct var- 
ious isothermal shock solutions, including twin shocks, two 
temperature shocks and so forth. 

Chevalier & Imamura (1982) performed a linear stabil- 
ity analysis for a one-dimensional shock driven by a constant 
velocity piston. For a radiative cooling function A oc p^T"', 
they found that radiative shocks were unstable for a ^ 0.4 
in a fundamental mode and unstable to overtone modes 
for a ^ 0.8. Bertschinger (1986) further extended their 
work to three dimensions and showed that transverse modes 
would be unstable for a ^ 1.0. Ryu & Vishniac (1987) ap- 
plied the linear theory to the dynamic instability of strong 
plane-parallel or spherical adiabatic blast shock waves in 
a gas medium with an initially uniform density and found 
that such blast waves to be unstable for an adiabatic in- 
dex 7 ^ 1.2. Blondin & Cioffi (1989) have shown that lo- 
cal instabilities are restricted to wavelengths less than the 
shock thickness. Vishniac (1993) examined the dynamical 
and gravitational instabilities of spherical shocks. Toth & 
Draine (1993) carried out linear stability analysis as well as 
numerical simulations to study effects of a transverse mag- 
netic fleld on magnetohydrodynamic (MHD) shock stability 
and found that a transverse magnetic fleld tends to stabilize 
the flow. MHD shocks in a gas medium of low fractional ion- 
ization are subject to a novel dynamical instability involv- 
ing the deformation of magnetic field lines (Wardle 1990). 
Decelerating radiative shock fronts (e.g., expanding SNRs 
in the snowplow phase) are subject to a ripping instability 
(e.g., Vishniac 1983; Vishniac & Ryu 1989). Drury (1984) 
noted that acoustic waves propagating towards a shock in 
the preshock medium may be amplified as they enter a re- 
gion with a steep cosmic-ray pressure gradient. Rayleigh- 
Taylor instabilities arise in the cosmic-ray precursor would 
further complicate the structure of the precursor when the 
density gradient due to these nonlinear sound waves opposes 
the cosmic- ray pressure gradient (e.g., Ryu 1993). 

Magnetic fields play important roles in various as- 
trophysical environments. Complex filamentary structures 
in molecular clouds, shapes and shaping of planetary 
nebulae, synchrotron radiation from SNRs, magnetized 
stellar winds, galactic winds, GRBs, dynamo effects in 
stars, galaxies, and galaxy clusters as well as other in- 
teresting problems all involve magnetic fields (e.g., Hart- 
mann 1998; Balick & Frank 2002). Due to the for- 
mation of filamentary structures in numerical simula - 
tions [e.g.. iPorter et al.) Jl994l) : iKlessen fc Burker^ (l200Cl) : 
lOstriker et alj i200ll) ] and the appearance of numerous fila- 
ments in observations (e.g., Falgarone et al. 2001), analyses 
on processes o f filament formation and evolution have also 
been pursued (K awach^^^lMiawa^98l: iHennebellg l2003l : 
iTillev fc Pudrit j2003l:IShadmehrill2005D . For these reasons, 
MHD shocks are under extensive investigations. Bazer fc Er- 
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icson (1959) were among the first to study hydromagnetic 
shocks for astrophysical applications. Whang (1984) stud- 
ied forward-reverse hydromagnetic shock pairs in the helio- 
sphere. Chakrabarti (1989, 1990) investigated MHD shocks 
in accretion discs. Numerical simulations were carried out 
to investigate the spherically symmetric shock interaction 
of pulsar wind nebula in a SNR (e.g., van der Swaluw et 
al. 2001). Recently, Del Zanna et al. (2004) re-examined 
this problem in an axisymmetric geometry. Ouyed & Pu- 
dritz (1993) studied obUque MHD shocks in disc winds from 
young stellar objects (YSOs) to explain broad, blueshifted 
forbidden emission lines observed in these sources. Duncan 
& Thompson (1992) proposed very strongly magnetized neu- 
tron star might be the origin of soft SGRs and anomalous 
X-ray pulsars (AXPs). Magnetic fields play an important 
role in the formation of the three-ring structure around su- 
pernova 1987A (e.g., Tanaka & Washimi 2002). Gaensler et 
al. (2005) found an expanding radio nebula produced by a 
giant flare from the magnetar SGR 1806-20 and interpreted 
it as ejecta colliding with pre-existing shells. 

In addition to Newtonian shocks, relativistic shocks are 
also investigated in various contexts. Blandford & McKee 
(1976) studied relativistic blast wave solutions in spherical 
geometry and calculate the active galactic nucleus (AGN) 
radiation spectrum using their blast wave solutions. Ken- 
nel & Coroniti (1984) considered ultra-relativistic MHD 
shock models for the Crab Nebula. Emmering & Chevalier 
(1987) studied relativistic MHD shocks and applied their re- 
sults to pulsar winds. Best & Sari (2000) re-examined the 
problem of Blandford & McKee (1977) and found second- 
type self-similar solutions to the problem of ultra-relativistic 
strong explosions. Perna & Vietri (2002) studied the prop- 
agation of a relativistic shock in an exponential atmo- 
sphere and applied their solutions to GRBs. Nakayama & 
Shigeyama (2005) reconsidered this problem in a plane- 
parallel geometry. Hidalgo & Mendoza (2005) considered 
imploding self-similar relativistic shock waves. Dynamics 
of relativistic magnetized blast wave was explored by Lyu- 
tikov (2002). Takahashi et al. (2002) investigated MHD adi- 
abatic shock accreting onto a rotating Kerr black hole. Sub- 
sequently, Fukumura (2004) studied isothermal shock for- 
mation around a Kerr black hole. Das et al. (2003) studied 
isothermal shocks around Schwarzschild black hole using the 
pseudo-Schwarzschild potential. Mildly magnetized internal 
relativistic shocks in GRBs have been invoked to explain the 
GRB prompt emission data (Fan et al. 2004). 

The role of a random magnetic field in our model anal- 
ysis is characterized by an important magnetic parameter A 
as defined by equation (|SJ . Different systems of astrophysi- 
cal objects involve a range of A values. The estimated A pa- 
rameter for the Crab Nebula is of the order of a magnitude 
around 10^ ~ 10®. For star forming regions, the estimated A 
is typically in the range of 0.01 ~ 0.1. For planetary nebulae, 
a typical A is approximately 10"'^. For a cluster of galax;ies, 
the estimated A is approximately 0.2. We would take A ~ 0.1 
typical for a star formation region for our later discussions. 
For a star formation cloud, we esitmate 



A~0.1 



Bn 



1.34 X lO-^G 



5 X 10-i9g/cm3 



2.24 X IQi'^cm 



This paper is structured as follows. The basic nonlin- 
ear MHD equations and the self-similar MHD transforma- 
tion are described in Section 2. In Section 3, we present the 
MHD shock conditions. In Section 4, we solve the nonlin- 
ear MHD ordinary differential equations (ODEs) numeri- 
cally to construct various MHD shock solutions. We provide 
comments and discussions in Section 5. Relevant technical 
details of mathematical analyses are summarized in several 
appendices for the convenience of reference. 



2 PHYSICAL ASSUMPTIONS AND 
BASIC MHD MODEL FORMULATION 

For a random tangled magnetic field on small scales, we for- 
mulate a large-scale MHD problem under the approximation 
of quasi-spherical symmetry. In spherical polar coordinates 
(r, 6, 0), the mass conservation equation is 



dp ^ 1 djr^pu) ^^ 
dt r^ dr 



(1) 



where u is the radial bulk fiow speed and p is the gas mass 
density. This mass continuity equation (0 above is equiva- 
lent to the following equations, namely 



dM 
dr 



4:Tvpr 



— — -a 
dt dr 



(2) 



where M{r, t) is the enclosed mass within r at time t. The 
isothermal radial MHD momentum equation is simply 



du du _ a^ dp GM 
dt dr p dr r^ 



Sirp dr " 



<Bl> 
Anpr 



(3) 



where a is the isothermal sound speed. By = {Bg + _B|)^ " 
stands for the random magnetic field component parallel to 
the MHD shock front and ~d^/dr = -GM{r,t)/r^ with 
^{r,t) being the gravitational potential. Here, < Bn > rep- 
resents a kind of ensemble average of random magnetic field. 
In our formulation, some terms involving cross correlations 
of radial and transverse magnetic field components have 
been ignored (see Appendix C of Yu & Lou 2005). Under 
the quasi-spherical symmetry, the Poisson equation relating 
the gas mass density p and the gravitational potential $ is 
automatically satisfied. Here in equation Q, we keep the 
magn etic tension force term that was ignored in equation 
(3) of IChiueh fc Choul Jl994D . 

By taking the electrical conductivity to be infinite (i.e., 
the ideal MHD approximation or complete frozen-in approx- 
imation), we arrive at the following relation 



B\\/{pr) ~ const 



(4) 



where the right-hand side is a constant of integration. The 
derivation of equation involves the mass conservation 
and the transverse components of the magnetic induction 
equation in the quasi-spherical approximation. For details, 
we refer the reader to Yu & Lou (2005). Using integral @, 
we define the dimensionless parameter A according to 



A = B?J{16nGp'^r^) 



(5) 
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which is a measure for the relative magnitudes of the mag- 
netic energy density and the self gravitational energy den- 
sity. In comparison with the previous work of Chiueh & 
Chou (1994), it is apparent that 



\ = (5{x)/{a'x') , 



(6) 



where a{x) and x axe the reduced density and the inde- 
pendent similarity variable defined immediately below in 
equation Q, respectively. Her e, variable I3(x) i s a di men- 
sionless function introduced bv lChiueh fc ChoxJ l|l994^ and 
stands for the reduced magnetic pressure, or equivalently, 
the reduced magnetic energy density. With this comparison, 
we know from their equation (11) that while our formula- 
tion differs from theirs in several aspects, the final magne- 
tosonic critical conditions are in fact identical. The ratio of 
the Alfven wave speed va to the isothermal sound speed a 
is given by 



VA 

a 



1/2 



with 



Bii 



(47r/9)V2 



The dimensionless independent similarity variable is defined 
by a; = r/{at) and the self-similar MHD transformation are 






■m{x) , u{r,t) = av{x) , 



ab(x) 

Vet 



(7) 



where the dimensionless a{x), m{x), v(x), <j){x) and b{x) are 
the reduced dependent variables for the mass density, the 
enclosed mass, the radial flow speed, the gravitational po- 
tential and the transverse magnetic field strength, respec- 
tively; and they are all dimensionless functions of x only. 
Note that h{x) = yXax and is related to the dimensionless 
function I3{x) of lChiueh &~C hou ( 1993) by b^ = p. 
It follows that equations Q reduce to 



[v — x)m + wi = 



(8) 



where the prime denotes the derivative with respect to x. 
Combining these two ODEs, we immediately obtain the fol- 
lowing two equations 

m = [x — v)x a , [x a{x — v)] = x a . (9) 

By equation (|5J, the physical requirement of a positive 
enclosed mass m{x) > is equivalent to the condition 
X — V > 0. Thus, solutions of v must lie to the upper-right 
of the straight line a; — w = in the plane —v{x) versus x. 

After substituting the self-similar variables into conti- 
nuity and momentum equations, we derive the following two 
coupled nonlinear MHD ODEs 

[(a; — n) — {1 + Xax )]v 

= {x-v) [a{x -v)- 2/x] , (10) 

\{x — v) — (1 + Xax )1 a /a 

= {x-v)[a-2{x-v)/x]+2Xxa . (11) 

Setting A = for the absence of the magnetic field, we read- 
ily recover the hydrodynamic formulation of Lou & Shen 
(2004). In the above two coupled nonlinear MHD ODEs JTUt 
and mi l, the magnetosonic critical curve is specified by 



X — V = 2/{ax) 



and 



x-v + 2Xx^ {x-vf . (12) 



Note that for A = and x — v > 0, equation 1121 becomes 
X — V = 1 and a — 2/x. Given the positiveness of a; — ii > 0, 
we come to the familiar condition x—v — 1 for the isothermal 
sonic critical line |Lou & Shcn 2004). Equation (a; — v)'^ = 



1 + Xax is equivalent to the equation (x 



= ia' + 



v\)/a'^, directly related to the fast magnetosonic condition 
as anticipated on the ground of physics. 

The asymptotic MHD solution behaviours of the cou- 
pled nonlinear MHD ODEs are summarized below. As MHD 
generalizations, these asymptotic solutions were derived in 
parallel to hydrodynamic results (see Appendix A of Lou & 
Shen 2004 for more details of their derivations) . We empha- 
size that there is no hydrodynamic counterpart for class HI 
asymptotic solutions as a; — > 0^ described below. In addition 
to v{x) and a{x), we also provide corresponding asymptotic 
solutions of P{x) for the reduced magnetic pressure for rele- 
vant physical information and a more complete presentation. 

In the limit of a: ^ +00, we have asymptotic solutions 



v{x) = V-f ■ 



A V 



+ — 
x^ 



(A/6 - l)(A - 2) + 2VV3 + A(2 - A)X/i 



(13) 



"(^) ^ -3 + 



A A{2~A) A{A-A)V 



2a;4 



+ 



3a"''' 



A{A - 3) (A/2 - 1) - (A - Q)AV'^/A + (2 - A)A^X/4: 



+ ••• , 
+ •■■ , 



(14) 



(15) 



where A and V are two constants of integration mainly for 
mass density and radial fiow speed, respectively. 

Class I: For a central MHD free-fall collapse in the limit 
of X ^ 0"^ , we have to the leading order 

V = -2F/x'-''^ - TT^^^^'^ In s - 2Lx'^''^ + ■■■ , (16) 



4F 



Q = F/x 



3/2 



-1/2 



SF 



In a; — Lx 



'1/2 



a \ 2 2 ^F , 
p — Xa X = h ■ 

X 



(17) 
(18) 



where F and L are two constants of integration. 

Class II: For the central LP-type MHD solutions being 
regular as a; -^ 0"*", we have 



2 (2 - 3I> - 18DA) 3 

^^^^D(2-3D-lsm^,^ 
18 

/3 = Xa X = \D X + ■ - • , 



(19) 

(20) 
(21) 



where D is an integration constant. 

Class HI: There exists a novel class of asymptotic MHD 
solutions which requires a sufficiently strong magnetic field 
and would disappear in the absence of magnetic field. As 
a: -^ 0"^, we write to the leading order 



v{x) — Hx + ■ 



(22) 



a(x) = K/x" + • 



(23) 



where K is an arbitrary integration constant and 

iJ= [2- A±(A^-4A)^''^]/2 < , (24) 

r; = (1 - _ff + 2A)/A with 2 < t? < 3 . (25) 



It is then straightforward to infer that as a: — > 0^ 



P{x) = Xa^x'^ 



k2('7-1) ' 



Bt 






(26) 



(27) 



For this new class of asymptotic MHD solutions, the mag- 
netic field strength scales as -B|| oc r~^''~^' at a fixed moment 
t, while at a given radius r, we have the scaling _B|| oc t''~^. 
The power index rj is such that —2 ^ —(j]^ 1) ^ ^1- By ex- 
pression 1241 , this type of MHD solutions exists only when 
the magnetic parameter A is greater than 4. 

It is possible for similarity MHD solutions to go across 
the magnetosonic critical curve smoothly. Along the magne- 
tosonic critical curve (e.g., Jordan & Smith 1977), the two 
MHD eigensolutions are governed by the following quadratic 
equation with z = ii', namely 



{v 



x\ 



xY 



-\- {x — v)z - 







(28) 



From equation 1281 , we obtain two types of MHD eigensolu- 
tions for z = v' along the magnetosonic critical curve. When 
the two roots of equation II28II are of opposite signs, type 1 
and type 2 eigensolutions are those with negative and posi- 
tive roots of equation 11281 , respectively. When the two roots 
of equation 1281 are of the same sign, type 1 and 2 eigen- 
solutions are those with smaller and larger absolute values 
respectively. In the open interval < a; < 2, type 1 and 
type 2 are exactly defined as such. When x > 2 (relevant to 
the LP-type MHD solution), dv/dx of type 1 has a larger 
absolute value, while dv/dx of type 2 has a smaller absolute 
value; that is, their magnitudes reverse for nodal points. In 
our current definition, no magnitude reversal would happen. 
When the point is a nodal point, type 1 and type 2 solu- 
tions remain always the smaller and larger ones respectively 
for the absolute value of dv/dx. To summarize, our defini- 
tion is not defined by the explicit expressions of dv/dx (e.g., 
1 — l/x* and 1/x,), but by their magnitudes and signs. It 
is easier to keep in mind their relevant physical properties. 



3 ISOTHERMAL MHD SHOCK CONDITIONS 
3.1 Perpendicular MHD Shock 

3.1.1 The One- Temperature Case 

The simplest type of MHD shock wave is the perpendicu- 
lar shock. In this case, the velocities of both the shock and 
plcisma are perpendicular to the magnetic field, which itself 
is parallel to the shock front. In a frame of reference mov- 
ing with the MHD shock front, the properties on both sides 
of the shock (pi, ui, _Bi,pi) and {p2,U2,B2,p2) are related 
by the equations for conservation of mass, momentum, en- 
ergy and magnetic fiux. In the isothermal approximation, 
we need not to consider the MHD energy equation. 
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In dimensional form, the perpendicular MHD shock 
conditions are (e.g.. Priest 1982) 



^21*2 = piui , 

P2 + -82/(871) + P2U2 

B2U2 — BlUi . 



- pi + Bf/{8n) + piul 



(29) 

(30) 

(31) 

In equation 13U1 . the term B^ /(Sn) represents the magnetic 
pressure. In equation 1311 . the quantity Bv gives the rate 
at which the magnetic fiux is transported across a unit sur- 
face area; physically, this quantity is actually proportional 
to the electric field component tangential to the shock front 
and thus should remain continuous across a shock. After 
straightforward manipulations, the solution to the above set 
of MHD shock equations can be expressed in terms of the 
mass density ratio P2/ p\ = A^, the upstream Mach num- 
ber M\ = u\/a, and the upstream plasma beta parameter'^ 
/3i = pi/{Bi/8tv) — la^ /v\^. The results are simply 

U2/UX = 1/X , 

B2/B1 = X , 

P2/P1 = X , 

where X is the positive root of the following quadratic equa- 
tion (see Appendix A for detailed derivations) 



/(A) = A" + (/3i + 1)A - l3iM? 







(32) 



this equation readily indicates the existence of only one pos- 
itive root of A. Physically, we should further require A > 1 
for a fast magnetosonic shock. 

When upstream conditions are specified, we can de- 
termine the downstream conditions behind a magnetosonic 
shock systematically. We simply take p2 as pd and pi as pu 
such that 



U2=Ud- 



U\ = Uu 



Ud = avd , 



A = P2/P1 = pd/pu , 

where Us is the radial shock speed, and subscripts u and 
d denote physical variables associated with upstream and 
downstream sides, respectively. It then follows that 

Qd/ou = A , 

{Vd - Xs)/{Vu - Xs) = 1/X , 
Ml =Vu ~ Xs , 

/?i = 2/{Xx^,a„) . 

Here, the value of A should be larger than 1 for a magne- 
tosonic shock - a kind of MHD fast shock. 

The analysis of quadratic equation 1321 shows that the 
two real roots of A have opposite signs. We naturally pick 
out the positive root among the two. 

When /(A = 1) < 0, i.e., 2 + f3i< (3iMf or 1 + 2/(3i < 
Ml, it has a positive root A > 1. This conclusion can also 
be reached by directly solving the quadratic equation and 

^ We emphasize here that /9i is the upstream plasma parame- 
ter, not the dimensionless I3{x) function for the reduced magnetic 
pressure defined in the previous context (see equation 151. 
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by requiring the positive root X > 1. Physically for a mag- 
netosonic shock to exist, the upstream flow speed relative 
to the shock speed must exceed the upstream fast magne- 
tosonic speed (a^ + Uai)'^ • This serves as a guide in con- 
structing magnetosonic shocks. 

When f{X = 1) > 0, i.e., 2 + /3i > f3iM^ or 1 + 2//3i > 
Ml, there is a positive root X < 1. 

Using the definitions of /3i and Mi in the condition 
1 + 2//3i = Ml and dropping the subscript 1, this equation 
describes exactly the behaviour of the magnetosonic critical 
curve {x — v)^ = 1 + Xax^. 



3.1.2 Two- Temperature MHD Shocks 

When the constant temperatures on the two sides of a mag- 
netosonic shock are different, we have the following two 
MHD jump conditions in the shock framework of reference 



ctd{vd — Xsd)ad — au{vu — Xsu)au , 
a^ [\aj_x^a/'2 + ad + ad{vd - x^d) ] 



(33) -^ 



(34) 



When the upstream conditions are given and we are going 
to determine the downstream conditions, we take p2 as pd 
and pi as p^. It then follows that 



ad/tSu = T , Xsd/xsu = l/r , 

dd/ciu = [Vu - Xsu)/[T{vd - Xsd)] = X , 

/3i = 2/{\xl^au) , 

X'^/Pi + T^X + Mi/X = l//3i + 1 + Mi 



(35) 
(36) 
(37) 
(38) 



Equation l|38^ is a cubic equation in terms of X, namely 

f{X) = X'-" + f3iT^X^ - {1 + I3i+ f3iM!)X + l3iM^ ^ .{39) 

For a physical magnetosonic shock solution, we should re- 
quire both r > 1 and X > 1. In constructing magnetosonic 
shock solutions, we then choose the positive root(s) X > 1 
among the three roots of cubic equation II38II for X. Since 
f{X = 1) = /3i(r2 - 1) > and f{X = 0) = f3iM? > 0, 
cubic equation H38^ must have at least one negative X root 
because /(— oo) = — oo and /(-foo) — -(-oo. The remaining 
two roots could be either both complex or both real. When 
the other two roots are either complex roots or both less 
than 1, there is no magnetosonic shock solution. When these 
other two real X roots are both larger than 1, there exist 
two possible magnetosonic shock solutions. From equation 
^, we derive f'{X) = df/dX as 



f'{X) = 3X^ + 2l3ir^X - (H- /?i -f l3iM^) 



(40) 



According to expression 1401 , the two real roots (one positive 
X+ and the other negative X-) of f'{X) — are given by 



X± = {-/3ir' ± [/3i\* + 3(1 + I3i+ l3iM^)]'^^}/3 



(41) 



By definition ^ for f{X) and the fact of f{X = 0) = 
PiMi > 0, it is clear that /(X_) should remain positive. For 
the situation of f{X+) > 0, the two remaining roots form a 
complex conjugate pair and there is no magnetosonic shock 
solution. For the situation of f{X+) ^ 0, the two remaining 
roots are real and positive. It is found that if the two roots 
are real only two situations can happen: (i) the two roots 



are both positive and less than 1; and (ii) the two roots are 
both larger than 1. The situation that one root is greater 
than 1 and the other positive root is less than 1 cannot 
occur because this would demand f{X = 1) < 0. 

We have just analyzed several solution properties of the 
cubic equation for the mass density ratio X given upstream 
physical conditions of a magnetosonic shock. However, in 
practical calculations, we did not search for root X > 1 
using the above procedure with r > 1. We adopted an al- 
ternative yet equivalent procedure described below. In fact, 
our later two-temperature calculations are performed with 
the following procedure involving one root of an equivalent 
cubic equation. 

Reciprocally, when downstream conditions are given, we 
can calculate the upstream conditions across a magnetosonic 
shock. Under this situation, we simply take p2 as pu and pi 
as pd and so forth. It then follows that 

" f, ^ = i, (42) 



ad 



ad 



Xsd 



{Vd ~ Xsd) 



= x 



Ml =Vd — Xsd , 

J3i = 2/{\xldad) , 

XV/9i + r^X + Mi/X = l//3i + 1 + Mi 



(43) 



(44) 



This is a cubic equation in X of exactly the same mathe- 
matical form as equation 13911 . namely 

f{X) = X-' + Pif''X^-{l + Pi+PiM?)X + PiM?=0.i'i5) 

For a physical magnetosonic shock, we naturally require 
f < 1 and X < 1. In numerical calculations, we choose 
the positive X root that is less than 1 among the three 
roots of cubic equation (I45I I. It turns out to be much sim- 
pler in analyzing this cubic equation. Because here /(X = 
1) = /3i(f2 - 1) < and /(X = 0) = /3iMf > 0, this cu- 
bic equation must then have three real roots, one negative 
and two positive. One positive root is larger than 1 and the 
other positive root is smaller than 1. In comparison with 
the cases of r > 1 in equation I39II . there is only one X root 
of cubic equation I45I I physically corresponding to a magne- 
tosonic shock solution. That the roots of the r > 1 cases in 
equation 1391 1 are more complex is related to the fact that 
/(X = 1) = Pi{t'^ ~ 1) > 0. Fortunately, in our specific 
model calculations, we just compute the f < 1 cases and 
the amount of our numerical computations can be greatly 
reduced. 

We have also considered the case of oblique MHD shocks 
in more details (not completely shown here) and find that 
oblique MHD shocks do not exist in our current model 
framework. The main reason that we cannot have oblique 
MHD shocks in our formalism is related to the frozen-in 
condition @. Thus from the two following shock relations 
for an obhque MHD shock (Priest 1982) 



B; 



2y 



{ul - vli)X 



Biy (u\~Xv\i) 

and 

X^Pl, 
Pi 

together with By oc p, we come to the conclusion that 
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{v 



i)X 



Xv\,) 



= x 



There is only a trivial solution with X = 1. From now on, we 
ignore the radial magnetic field component and just consider 
perpendicular MHD shocks. 




4 SELF-SIMILAR ISOTHERMAL 
MHD SHOCK SOLUTIONS 

In this section, we present and discuss properties of semi- 
complete global solutions with perpendicular MHD shocks. 



4.1 Extensions of Previous Shock Results by 
Including a Random Magnetic Field 

With the inclusion of a random magnetic field, we can 
extend previous hydrodynamic similarity shock solutions. 
More specifically, we would extend the results of Tsai & Hsu 
(1995), Shu et al. (2002) and Shen & Lou (2004). For a ran- 
dom magnetic field, we envision a simple "ball of thread" 
scenario in a huge spatial volume of gas medium. A mag- 
netic field line follows the 'thread' meandering within a thin 
spherical 'layer' in space in a random manner. In the strict 
sense, there is always a random weak radial magnetic field 
component such that random magnetic field lines in adja- 
cent 'layers' are actually connected throughout in space. By 
taking a large-scale ensemble average of such a magnetized 
system, we are then left with 'layers' of random magnetic 
field components transverse to the radial direction. Very re- 
cently, Bian & Lou (2005) carried out extensive investiga- 
tions on self-similar isothermal shock solutions in the ab- 
sence of magnetic field, while under the approximation of 
a quasi-spherical symmetry in the sense described above, 
Yu & Lou (2005) constructed smooth self-similar isothermal 
MHD solutions involving a random magnetic field but with- 
out shocks. To a greater extent, our present investigation is 
to combine the analyses of Bian & Lou (2005) and Yu & 
Lou (2005) as well as to explore new possibilities associated 
with magnetic field. 



4.1.1 MHD Extensions for Shocks of Tsai & Hsu 

Alternative to the results of Larson (1969a, b) and Pension 
(1969a, b), Shu (1977) constructed the EWCS and devel- 
oped an inside-out collapse scenario for the process of form- 
ing low-mass protostars (e.g., Shu et al. 1987). Tsai & Hsu 
(1995) considered a self-similar shock travelling into a static 
SIS envelope characterized by r = 1, Xsd = a;su = Xs, Uu = 
and Qu = 2/x^. Their shock connection condition becomes 
Vd = Xs — 1/xs and a^ ~ 2 (see their equation 16). Their 
global shock solutions have different asymptotic behaviours 
near the origin (see their fig. 6). The diverging free-fall solu- 
tion near the origin is the Class I solution. The converging 
LP- type solution near the origin is the Class II solution. In 
our formalism, we can generalize the shock results of Tsai 
& Hsu (1995) by incorporating a random magnetic field. 
When a magnetic field is included, the MHD Class I shock 
solution generalizing that of Tsai & Hsu (1995) crosses the 
magnetosonic critical curve at x, — 0.0369 with the shock 
location at Xs = 1.3599 for an outgoing MHD shock front 
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_^^^^^^^MTH 






MEWCS / ~~~" 
MLP 


D=11.26 


_i 



^^"""^^^^ir 



Figure 1. The MHD case of A = 0.1. The MHD extensions of 
Tsai & Hsu (MTH) type similarity shock solutions (Class I MHD 
shock with free-fall diverging behaviour near the origin and Class 
II MHD shock with a finite reduced density near the origin and 
an eigenvalue D = 11.26). The negative reduced radial velocity 
—v{x) (upper panel) and the reduced density a{x) in logarithmic 
scale (lower panel) versus x in linear scale are shown respectively. 
MHD Expansion-Wave Collapse Solution (MEWCS) as well as 
MHD LP (MLP) solution are also displayed. Dash-dotted line 
represents the magnetosonic critical curve. Class I MHD Tsai & 
Hsu (MTH) type similarity shock location is at Xs = 1.3599 and 
this solution passes through the magnetosonic critical curve at 
x, = 0.0369. Class II MHD Tsai & Hsu (MTH) type similarity 
shock location is at Xs = 1.4194. 



travelling at a constant speed of 1.3599 times the isother- 
mal sound speed a. The corresponding central mass accre- 
tion rate is 0.0726. For the MHD Class II shock solution, a 
central expansion with a finite central density (LP-type so- 
lution) matches with a magnetostatic SIS envelope across an 
MHD shock. The shock is located at Xs = 1.4194 and thus 
shows a higher shock speed of 1.4194 times the isothermal 
sound speed a and a reduced core density D = 11.26 [see 
asymptotic solutions 11911 . 1201 . and 1211 1. In Figure Q and 
Figure 121 we plot the MHD extensions of Tsai & Hsu type 
self-similar shocks. Note that these MHD generalizations of 
Tsai & Hsu shock solutions are just two special cases among 
magnetized similarity shock solutions into a magnetostatic 
SIS envelope. Details of the procedure for constructing these 
MHD shock solutions will be discussed in section 4.2. 



4.1.2 MHD Extensions for Shock Solutions of Shu et al. 

Shu et al. (2002) extended the Class II shock solution of 
Tsai & Hsu (1995) and used these solutions to model 'cham- 
pagne fiows' of HII regions surrounding massive OB stars. 
In a similar manner, by varying the D and Xs parameters 
with V = and a fixed A, MHD extensions for shock solu- 
tions of Shu et al. (2002) are also constructed and presented 
in Figure 13 and Figure 2] Here the parameters D, V, \ and 
A, K and H in later discussion are all integration constants 
in asymptotic solutions 1131 to 1231 . Both the shock speed 
and shock strength increase with a decreasing D parameter. 
In the limit of Z> — » 0"*" numerically, the mass parameter 
A approaches 0^ accordingly and the reduced radial speed 
v{x) converges to an invariant form (Shu et al. 2002; Shen 
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Figure 2. The MHD case of A = 0.1. The MHD extensions of 
Tsai & Hsu (MTH) type similarity shock solutions. The ratio of 
Alfven wave speed to the isothermal sound speed va/ci versus x is 
shown in linear scales. MHD Expansion- Wave Collapse Solutions 
(MEWCSs) and MHD LP (MLP) shock solutions are displayed. 



& Lou 2004) with an invariant fastest and strongest MHD 
shock located at x^ = 2.56. In other words, as D — > 0^ nu- 
merically, behaviours of the MHD shock solutions gradually 
become the same as those of the hydrodynamic shock solu- 
tions, and for the same parameter D, the MHD shock speeds 
are faster than the hydrodynamic shock speeds. These shock 
solutions are MHD generalizations of 'champagne breezes' 
(see curves in our Fig. El and the heavy solid curve in our 
Fig. 13 and compare with fig. 1 of Shu et al. 2002). In Table 
1, we summarize properties of the displayed shock solutions 
in the semi-complete space < a; < -l-oo. 

Details of the procedure for finding these semi-complete 
self-similar MHD shock solutions will be discussed presently 
in section 4.3. 



^.1.3 MHD Extensions for Shock Solutions of Shen & Lou 

Shen & Lou (2004) further extended the Class I shock solu- 
tions by matching with various asymptotic fiows and mod- 
elled the dynamical evolution of young stellar objects such as 
Bok globule B335 system to account for the observationally 
inferred mass density and fiow speed profiles as well as the 
estimated central mass accretion rate. When a random mag- 
netic field is included, MHD extensions for this type of shock 
solution can also be constructed. The downstream side of 
such an MHD shock is part of the first type 2-type 1 solution, 
i.e., the MHD EECC solution (MEECC; Yu & Lou 2005) 
which crosses the magnetosonic critical curve analytically. 
The two magnetosonic critical points are at a;*(l) = 0.103 
and a;, (2) — 1.811, respectively. By choosing different shock 
locations Xs = 0.4, 0.9, 1.5018, 1.57 as in Figure|Sland Fig- 
ure|S| we readily construct various upstream MHD solutions 
as a; ^ +oo. Magnetized 'champagne flows' of HII regions 
around massive OB stars can also be constructed by allowing 
for MHD flows at large x. Specific examples of such solutions 
are presented in Figure [7| and Figure |H| Here, the upstream 
flows are not necessarily limited to a static magnetized SIS 
envelope or magnetized 'champagne breezes' with V^ = 0. 
With different MHD shock locations or speeds. Class II sim- 




Figure 3. The MHD generalization for similarity shock solutions 
of Shu et al. (2002) with A = 0.1. For self-similar MHD shock so- 
lution curves from bottom to top, D = 2/3, 1, 1.67, 4, 11.26, 
respectively. The negative reduced radial speed —v{x) in linear 
scale (upper panel) and the reduced mass density a{x) in loga- 
rithmic scale (lower panel) versus x in linear scale. Both the MHD 
shock speed and strength increase with a decreasing D parame- 
ter. As D approaches 0+, the reduced speed v{x) converges to 
an invariant form with an invariant fastest and strongest shock 
located at Xs = 2.56 (see also the heavy solid curve in Fig. 171. 
Dashed line represents the magnetosonic critical curve. 




0.5 1 



Figure 4. MHD extensions for similarity shock solutions of Shu 
et al. (2002) with A = 0.1. The ratio of Alfven wave speed and 
isothermal sound speed va/ci versus x is shown. Both the shock 
speed and strength increase with a decreasing D parameter. 



ilarity MHD solutions can be matched with either asymp- 
totic MHD outflows {V > 0) or asymptotic MHD inflows 
{V < 0). 



4.2 MHD Similarity Shocks into 
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Table 1. MHD Extensions for Shocks of Shu et ah (2002) with A = 0.1. 



description 


A 


Xs 


Vd 


"d 


Vu 




Ctu 


D = 11.26 


2.00 


1.42 


0.532 


1.59 







0.992 


_D = 4 


1.74 


1.69 


0.837 


1.14 


0.230 




0.666 


D = 1.67 


1.42 


1.93 


1.11 


0.816 


0.423 




0.441 


D = 1 


1.20 


2.06 


1.27 


0.653 


0.534 




0.337 


D = 2/3 


1.02 


2.15 


1.38 


0.537 


0.620 




0.269 


D = 10-5 


4.30 X 10-5 


2.56 


1.91 


2.10 X 10-5 


1.02 


8 


80 X 10-" 


D = 10-^ 


4.30 X lO-fS 


2.56 


1.91 


2.10 X 10-6 


1.02 


8 


80 X lO-'^ 



These solutions can be viewed as Class II similarity MHD shock solutions matched with 
asymptotic MHD breezes. Columns 1 to 7 provide relevant parameters for MHD shock 
solutions: D is the key parameter in the LP- type asymptotic MHD solutions I19i . 1201 
and I21i at small x; A is the upstream mass density parameter; Xs indicates either MHD 
shock location or MHD shock speed; v^ is the reduced speed downstream of an MHD 
shock; aj is the reduced density downstream of an MHD shock; Vu is the reduced speed 
upstream of an MHD shock; and au is the reduced density upstream of an MHD shock. 




x.=0.103 




Figure 5. The first MHD EECC (MEECC) solution with A = 0.1 
crossing the magnetosonic critical curve analytically at x«(l) = 
0.103 and a;* (2) = 1.811 [within a;* (2) is light solid curve, while 
outside Xt,(2) is dotted curve]. The light solid curves represent 
Class I MHD shock solutions with the downstream as part of 
the first MEECC solution and with shock location Xs at 0.4, 0.9, 
1.5018, 1.57, respectively; the heavy solid curves represent Class 
II MHD shock solutions with the reduced mass density of the 
central core D = 0.1 and with the MHD shock location Xs at 1.30, 
2.00, 2.80, respectively. Dash-dotted line represents magnetosonic 
critical curve. Dashed line is the mass demarcation line x — v = 0. 



a Magnetostatic SIS Envelope 

MHD extensions for shock solutions of Tsai & Hsu (1995) 
shown earlier are just two cases of our following similarity 
MHD shock solutions into a magnetostatic SIS envelope. 
In Tsai & Hsu (1995) and Shu et al. (2002), two classes 
of shock solutions were constructed to match with a static 
SIS envelope with V = and A = 2 or A < 2. In Bian 
& Lou (2005), hydrodynamic shock solutions of Class I and 
Class II were extensively explored to match with a static 
SIS envelope. One can also follow the procedure of Hunter 
to derive the Class I and Class II MHD shock solutions in a 
parallel manner. Here, we would systematically explore the 
Class I and Class II MHD shock solutions matched with a 
magnetostatic SIS envelope by surveying possible solutions 



Figure 6. The ratio of the Alfven wave speed to the isothermal 
sound speed va/o, versus x with A = 0.1 corresponding to the case 
of Fig. |5] The light solid curves represent Class I MHD shock so- 
lutions with the downstream as part of the first MEECC solution 
and with the shock location Xs at 0.4, 0.9, 1.5018, 1.57, respec- 
tively; the heavy solid curves represent Class II MHD shock solu- 
tions with the reduced mass density of the central core D = 0.1 
and with the MHD shock location Xs at 1.30, 2.00, 2.80, respec- 
tively. 



in the speed-density phase diagram of v and a. In Figure 
1101 we present the Class I MHD shock solutions. In Figure 
1141 we provide the Class II MHD shock solutions. In Table 
2, we summarize relevant parameters of the first three MHD 
shock solutions of both Class I and Class II. 

We first consider Class I MHD shock solutions with an 
outer magnetostatic SIS envelope, having downstream solu- 
tions diverging towards the centre and upstream solutions 
of mass parameter A = 2. Tsai & Hsu (1995) obtained a 
hydrodynamic shock solution for such a case. Bian & Lou 
(2005) systematically extended the shock results of Tsai & 
Hsu (1995) and obtained a wide variety of similarity shock 
solutions. We shall fully explore the speed-density phase di- 
agram for constructing semi-complete MHD shock solutions. 
For every assigned value for x«(l) [a;*(l) < Xm. where Xm is 
a chosen meeting point] and after integrating towards the 
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Figure 7. The case of A = 0.1. The negative reduced radial 
speed —v(x) (upper panel) and the scaled reduced mass denstiy 
R(x) = x^oi.{x)/ D (lower panel) versus x. The solid curves form 
the family of Class II MHD shock solutions in the invariant form 
with the reduced mass density of the central core D — > 0+; the 
heavy solid curve is the MHD shock 'champagne breeze' solution. 
Dash-dotted curve in the upper panel is the magnetosonic critical 
curve. 




Figure 9. The phase diagram of Class I similarity MHD shock 
solutions with A = 0.1. The meeting point is chosen at Xm = 
0.5. The type 2 MHD eigensolution is used at the point x*(l) to 
integrate towards the meeting point. In the range oi x > Xs, the 
magnetostatic SIS envelope with A = 2 and V = is used to 
get this diagram. The values of Xt{l) and Xs of the first three 
intersection points are [xt(l) = 0.0369, Xs = 1.3599], [a;t(l) = 



2.044 X 10" 



1.1266], respectively. 



Xs = 0.7807] and [x*(l) = 1.10 X 10" 
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Figure 8. The ratio of the Alfven wave speed to the isothermal 
sound speed Vj\/a versus x with A = 0.1, corresponding to MHD 
shock solutions in Fig. 171 



meeting point Xm. with a type 2 MHD eigensolution, we ob- 
tain a pair of {v, a\ at the meeting point Xm in the so-called 
phase diagram of v versus a. For a sequence of such a;.(l) 
values, a series of {v, a} pairs is obtained, giving rise to 
a curve in the phase diagram. The upstream is part of a 
MEWCS. By choosing an MHD shock location Xs > Xm, we 
determine the corresponding upstream values Vu and Qu at 
Xs along the MEWCS solution. We then calculate Vd and 
Qd using the isothermal MHD shock conditions in terms of 
the upstream values v^ and «„. Starting from Vd and ad 
at Xs, we then integrate the coupled nonlinear MHD ODEs 
backwards to the meeting point Xm to obtain another pair 
of {v, q}. In other words, every value of Xs corresponds to a 
pair of {v, ct}. For a sequence of Xs values, another curve in 
the phase diagram is thus produced. The intersection points 



of the two phase curves represent matches of this type of 
similarity MHD shock solutions. 

In Fig. 1^1 we show the relevant phase diagram of v and 
a following the matching procedure described above for a 
chosen meeting point at Xm = 0.5. Relevant parameters of 
the first three intersection points in the phase diagram are 
(x,(l),2;,)=(0.0369, 1.3599), (2.044 x lO^**, 0.7807), (l.lOx 
10"'*, 1.1266), respectively. In Figs. llUl and llll we display the 
first three global MHD shock solutions of this type. Enlarged 
portions of these solutions are given in Fig 1121 

When it comes to the Class II MHD shock solutions 
that connect downstream MHD LP solutions with a magne- 
tostatic SIS envelope, the method of finding the solutions is 
similar to those described above. The key difference is that 
the downstream is now replaced by the MHD LP type solu- 
tions. In parallel, we could obtain similar phase diagram to 
identify the intersection points of phase curves. In Fig. 1131 
we present a relevant phase diagram of v and a to match 
similarity MHD shock solutions. In Fig. 1141 and Fig. 1151 we 
display the first three global MHD shock solutions of this 
type. 

In Table 2, we summarize the properties of Class I and 
Class II MHD shock solutions matched with a magnetostatic 
SIS envelope, i.e., MEWCS solutions. 



4.3 Similarity MHD Shock Breezes 

MHD extensions for shock solutions of Shu et al. (2002) de- 
scribed in previous section should be viewed as a subset of 
Class II type MHD shock solutions. We now consider Class I 
type MHD shock solutions and allow the mass parameter A 
to vary (i.e., A cannot be equal to 2 in general) for a given x* 
corresponding to a specific reduced central mass accretion 
rate mo. Both parameters A and Xs are adjusted gradually 
to match the upstream and downstream solutions in the 
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Table 2. Semi-complete Class I and II similarity MHD shock solutions matched with a 
magnetostatic SIS envelope. 



Description 


mo 


Xs 


Vd 


ad 


Vu 


OLu 


N 


a;.(l) = 0.0369 


0.0726 


1.3599 


0.4421 


1.6024 





1.0815 


1 


a;*(l) = 2.044 X 10"'' 


4.09 X 10-4 


0.7807 


-0.2450 


2.6333 


-0.3428 


2.4019 


2 


x*(l) = 1.10 X 10-** 


2.20 X 10-** 


1.1266 


0.0567 


1.6593 





1.5758 


3 


D = 11.26 


- 


1.41944 


0.5318 


1.5847 





0.9926 





D = 2156 


- 


0.77598 


-0.2478 


2.6593 


-0.3501 


2.4177 


1 


D = 4.04 X lO'^ 


- 


1.12668 


0.0569 


1.6593 





1.5755 


2 



Columns 1 to 8 provide relevant parameters for MHD shock solutions: a;,(l) is the inner 
magnetosonic critical point for Class I MHD solutions and D is the central 'density 
parameter' of MHD LP type symptotic solution for constructing Class II MHD solutions; 
mo is the reduced central mass accretion rate; Xs corresponds to both the shock location or 
shock speed; v^ is the downstream reduced speed; a^ is the downstream reduced density; 
Vu is the upstream reduced speed; oid is the upstream reduced density and the number 
of stagnation points A'^ indicates the number of nodes for sub-magnetosonic self-similar 
oscillations in MHD shock solutions. 



;.=0.1, A=2, V=0 

1 St solution : x.(1 )=0.0369, x =1 .3599 
2nd solution : x.(1)=2.044e-4, x =0.7807 
3rd solution : x.(1)=1.10e-6, x =1.1266 
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Figure 10. The first three Class I MHD shock solutions with 
A = 0.1. The MHD shock solution 2 is an accretion MHD shock. 
The reduced radial speeds of the preshock and postshock are both 
negative. This similarity accretion MHD shock expands at a con- 
stant sub-magnetosonic speed. The type 2 MHD eigensolution is 
chosen at the point Xt(l) to integrate towards the meeting point. 
In the range of a; > Xs, the magnetostatic SIS envelope with 
A = 2 and y = is used to get this diagram. The dash-dotted 
curve represents the magnetosonic critical curve. The dashed line 
\s x — v = 0. The corresponding speed ratio v^/a is shown in Fig. 

ini 



phase diagram of v versus a. We integrate from a specified 
Xt along the magnetosonic curve (in our calculations, we 
choose X, = 0.0369, 0.103, 0.2, 0.4, 0.7 as several examples 
of illustration) with a type 2 MHD eigensolution towards the 
shock front at Xs to obtain the two parameters Vd and ad- 
By the isothermal MHD shock conditions, we determine «„ 
and Qi, for a given pair of Vd and ctd- Starting from Vu and 
a-u at Xs, we then integrate further the coupled nonlinear 
MHD ODEs to a chosen meeting point Xm = 3.0. Here, we 
apply the analytical asymptotic MHD solution (i) as given 
by equations (I13II — I15II at a; = 20 as the far-away 'bound- 
ary condition' to integrate backwards to Xm ~ 3.0. For a 
given a;* value and two series of Xs and A values, we then 



Figure 11. Corresponding to the shock solutions in Fig. 1101 the 
ratio of the Alfven wave speed to the isothermal sound speed va/ci 
versus x with A = 0.1. The first three Class I MHD shock solutions 
identified in the a — v phase diagram of Fig. El The type 2 MHD 
eigensolution is taken at the point x«(l) to integrate towards the 
meeting point. In the range of a; > Xs, the magnetostatic SIS 
envelope with A = 2 and V = is used for these MHD shocks. 



draw a relevant phase diagram to determine the intersec- 
tion point of the two phase curves, i.e., Xs and A. We could 
then construct such a similarity MHD shock breeze using 
the values of Xs and A at the intersection point. In Fig. 1161 
we display a sample phase diagram for the case of a;, = 0.2. 
For a;* — 0.0369, 0.103, 0.4, 0.7 respectively, we can ob- 
tain qualitatively similar phase diagrams (not shown here to 
avoid redundancy) and their corresponding a;^ and A values. 
In Fig. 1171 we present Class I isothermal MHD shock breeze 
solutions. The corresponding ratios of the Alfven wave speed 
VA to the isothermal sound speed a are shown in Fig llSI Ta- 
ble 3 contains the major relevant parameters for the shock 
solution examples of illustration. 

Methods of determining Class II solutions of MHD 
shock breeze are similar to the above procedure. The main 
difference is to adopt the MHD LP type solutions for the 
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Figure 12. Enlarged portions of Class I MHD shock solution 
curves near x — » 0"*" of Fig. IIUI emphasizing the diverging and 
oscillatory behaviours of this family of MHD shock solutions as 
a; — > 0"*" . The a;— axis is shown in the logarithmic scale. The dashed 
line on the right is the magnetosonic critical curve. The undula- 
tory profiles of the curves marked by numerals 1, 2, 3 represent 
self-similar sub-magnetosonic radial oscillations with one, two, 
three stagnation points, respectively. 



>.=0.1, V=0, A=2 

1 St solution : X =1 .41 944, D=1 1 .26 
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2nd solution : x =0.77598, D=2156 
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3rd solution :x =1.12668, D=4.04x10° 




Figure 14. The first three class II MHD shock solutions with 
A = 0.1 identified in the phase diagram of Fig. lliji The shock solu- 
tion 2 is an accretion MHD shock with the reduced radial speeds 
of the preshock and postshock being both negative. This accretion 
MHD shock front expands at a constant sub-magnetosonic speed. 
The LP type MHD solution is used at near the origin to integrate 
outward to the meeting point. In the range oi x > Xs, the magne- 
tostatic SIS envelope with A = 2 and V = is specified for this 
diagram. The dash-dotted line represents the magnetosonic criti- 
cal curve and the straight dashed line represents the demarcation 
line of x — D = for a positive enclosed mass. The speed ratio of 
the Alfven wave speed to the isothermal sound speed va/o, versus 
X is shown in Fig. 1151 
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Figure 13. The phase diagram of Class II similarity MHD shock 
solutions with A = 0.1. The meeting point is chosen at Xm=0.5. 
The MHD LP type solution is used at the point near origin to in- 
tegrate outward to the meeting point. In the range oi x > Xs, the 
magnetostatic SIS envelope with A = 2 and V = is specified for 
this diagram. The values of Xs and D of the first three intersection 
points are (xs = 1.41944, D = 11.26), (xs = 0.77598, D = 2156) 
and (xs = 1.12668, D = 4.04 x 10^^), respectively. 
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Figure 15. For the first three Class II MHD shock solutions with 
A = 0.1 identified in the phase diagram of Fig. Iliji and shown in 
Fig. 1141 the corresponding ratio of the Alfven wave speed to the 
isothermal sound speed Vj\/a versus x is in display here. The 
MHD LP type solution is used near the origin to integrate out- 
ward to the meeting point. In the range of x > Xs, the magne- 
tostatic SIS envelope with A = 2 and V = is specified for this 
diagram. 



downstream. Results of Class II MHD shock solutions have 
been discussed already in subsection 4.1.2 and we shall not 
repeat here. 



4.4 Self-Similar Twin MHD Shock Solutions 

In parallel with the investigation of Bian & Lou (2005), we 
can also construct twin shock solutions in the presence of 
a random magnetic field. All MHD shock solutions shown 
so far contain just a single MHD shock. In this subsection. 
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Table 3. Semi-complete Class I similarity MHD shock solutions matched with an asymp- 
totic MHD breeze. 



Xt 




mo 


A 


Xs 


Vd 


Oid 


Vu 


a-u 


X* 


= 0.0369 


0.0724 


2 


1.3599 


0.4414 


1.6020 





1.0828 


Xt, 


= 0.103 


0.1952 


1.8530 


1.5018 


0.5805 


1.2907 


0.1507 


0.8802 


Xt, 


= 0.2 


0.3586 


1.8064 


1.5146 


0.5718 


1.2050 


0.2020 


0.8655 


Xt 


= 0.4 


0.6356 


1.8514 


1.4216 


0.4265 


1.2520 


0.1870 


1.0091 


Xt 


= 0.7 


0.9107 


1.9484 


1.2492 


0.1897 


1.4401 


0.1014 


1.3294 



Columns 1 to 8 provide the relevant parameters for MHD shock solutions: a;» is the 
magnetosonic critical point for Class I MHD shock breeze solutions; mg corresponds to 
the reduced central mass accretion rate; A is the mass density parameter of the upstream; 
the shock location is at Xs\ Vd is the downstream reduced velocity; ad is the downstream 
reduced density; Vu is the upstream reduced velocity; and au is the upstream reduced 
density. 



£5 0-204 



meeting point x =3, x.=0.2 

1.87* 




A=1.74 



> 

I 




Figure 16. An example of illustration for the phase diagram 
of Class I MHD shock breeze solution with A = 0.1 for the case 
of Xt = 0.2. For other values of x«, similar phase diagrams are 
necessary to identify the relevant global MHD shock solutions. 
The type 2 MHD eigensolution is specified at the point a:«(l) to 
integrate towards the meeting point. Analytical asymptotic MHD 
solution (i) as given by equations I13i — 1151 with V = and a 
varying A parameter is specified to integrate also towards the 
meeting point. 



Figure 17. The Class I MHD shock solutions with A = 0.1. The 
type 2 MHD eigensolution is specified at the point x«(l) along 
the magnetosonic critical curve to integrate towards the meeting 
point Xm. Analytical asymptotic MHD solution (i) as given by 
equations I13i — 1151 with V = and a varying A parameter 
is invoked to also integrate towards the meeting point Xm- The 
dash-dotted line represents the magnetosonic critical curve. The 
speed ratio of vj^/a is shown in Fig. 1181 



we are going to explore examples of semi-complete simi- 
larity solutions with twin MHD shocks. The upstream be- 
tween the outer MHD shock and the magnetostatic state 
are tangent to the magnetosonic critical curve at the point 
of x,{2) = (1 + 2\Y''^, which actually is part of MEWCS 
(Yu & Lou 2005). We may shift the location x, (2) somewhat 
towards the origin into the region of < a; < (1 + 2A)^" and 
the upstream solution of the inner MHD shock may cross the 
magnetosonic critical curve analytically with a type 2 MHD 
eigensolution. We can subsequently construct Class I and 
Class H twin MHD shock solutions in the a — v phase dia- 
gram. For constructing Class I twin MHD shock solutions, 
we first specify the value of a;* (2) and adjust the value of 
a;*(l) point and the shock location Xs(l) between a;*(l) and 
a;* (2) to search for similarity MHD shock solution crossing 
the magnetosonic critical curve analytically using the inter- 
sections of phase curves in the a — v phase diagram. More 
specifically, we choose a meeting point Xm = 0.5 and inte- 
grate an upstream solution of the inner shock from a;, (2) 



for Class I solutions using a type 2 MHD eigensolution as 
initial conditions to a chosen shock location Xa{V) in order 
to get the upstream values of «„ and q„ at the inner shock 
location Xs(l). We next use the isothermal MHD shock con- 
ditions to determine the downstream values of Vd and Ud at 
the inner shock location Xs{l). We then integrate from Vd 
and ad at a;s(l) towards the meeting point Xm ~ 0.5. Mean- 
while, we also integrate forward from a;«(l) again using type 
2 MHD eigensolution to Xm ~ 0.5. In Fig. 1191 we show the 
phase diagram of Class I MHD shock solution for two cases 
a::*(2) = 0.9 and 1.0, respectively (n.b., phase curves of Class 
H MHD shock solution are also included in this figure) . From 
the intersections of phase curves, we obtain the relevant pa- 
rameter pair [a;«(l), a;a(l)] of Class I MHD shock solutions 
as (1.44 X lO""*, 0.8352) for the case of a;.(2) = 0.9 and as 
(1.88 X 10"^, 0.7901) for the case of a;. (2) = 1.0. We show 
the corresponding Class I MHD shock solutions for —v ver- 
sus X in Fig. 1201 as an example of illustration. In Fig. 1211 we 
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Figure 18. The ratios of the Alfven wave speed and the isother- 
mal sound speed va/o, with A = 0.1 corresponding to the MHD 
shock breeze results in Fig. 1171 



show the corresponding ratios of the Alfven wave speed va 
to the isothermal sound speed a. 

For constructing Class II twin MHD shock solutions, 
we first specify the value of a;, and adjust the value of cen- 
tral reduced mass density D parameter together with the 
MHD shock location Xgil) between and x^ to search for 
similarity MHD shock solutions using the intersections of 
phase curves in the a — v phase diagram. From the intersec- 
tions of phase curves in Fig. 1191 we identify the relevant pa- 
rameter pair [D, 2:3(1)] of Class II MHD shock solutions as 
(3059, 0.8297) for the case of a:. = 0.9 and as (2344, 0.7854) 
for the case of s, = 1.0. We present the corresponding Class 
II MHD shock solutions for —v versus x in Fig. 1221 as an- 
other example of illustration. The corresponding ratios of 
the Alfven wave speed va to the isothermal sound speed a 
are shown in Fig 1231 

In Fig. 1191 relevant phase curves of both Class I and 
Class II twin MHD shock solutions are presented, where the 
meeting point Xm ~ 0.5 is chosen. In this phase diagram, 
both cases of a;, (2) — 0.9 and x„{2) = 1.0 for Class I solution 
and of x^ = 0.9 and a;, = 1.0 for Class II solution are 
plotted. Using the intersection points of phase curves in this 
phase diagram, we can readily construct the inner Class I 
and Class II MHD shock solutions, respectively. 

The above procedures are implemented to search for 
the inner MHD shock solutions. In the range of a; > x,(2) 
for Class I solution (or a; > a;* for Class II solution), we 
can further construct the outer MHD shock by choosing 
different outer shock location a;s(2). For a;, (2) — 0.9 (for 
Class I) or x^, — 0.9 (for Class II), we choose a;s(2) = 
0.95, 1.00, 1.10, 1.15, 1.20, 1.25, respectively, while for 
a;* (2) = 1.0 (for Class I) or x, — 1.0 (for Class II), we choose 
a;s(2) — 1.05, 1.10, 1.15, respectively. The upstream asymp- 
totic MHD solutions across the outer shock can match with 
the analytical asymptotic solution I13II — 1151 1 for x -^ +00. 
When the value of the outer shock location Xs{2) is smaller, 
the upstream would be an MHD inflow (contraction or accre- 
tion). As the value of MHD shock location a;s(2) increases, 
the upstream would become an MHD outflow (expansion 
or wind or breeze). Relevant solutions are displayed in Fig. 

Hl-Fig.llTl 



Figure 19. A phase diagram of v versus a at a chosen meeting 
point Xm = 0.5 for constructing MHD shock solutions with A = 
0.1. In this phase diagram, both cases of a:, (2) = 0.9 and a;* (2) = 
1.0 for Class I MHD shock solutions and of a;» = 0.9 and a;» = 1.0 
for Class II MHD shock solutions are shown. 
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Figure 20. Class I MHD shock solutions with A = 0.1 for 
Xt(2) = 0.9 and 1.0 shown in terms of —v{x) versus x by the 
heavy and light solid curves, respectively. The two correspond- 
ing values of inner magnetosonic point a't(l) are 1.44 X 10~* 
and 1.88 X 10~^, respectively. The two corresponding inner MHD 
shocks are located at Xs = 0.8352 and Xs = 0.7901, respectively. 
The dash-dotted line represents the magnetosonic critical curve. 
The speed ratio va/o, versus x is shown in Fig. 1211 



4.5 Novel Class III MHD Shock Solutions 

In our numerical investigation, we also construct a unique 
Class HI solutions exclusive to MHD shock solutions. This 
type of solutions carry a mixed feature of both Class I and 
Class II solutions. Their reduced velocity is a finite conver- 
gent inflow, similar to Class II solutions towards the cen- 
tre. Their reduced mass density is divergent, characteris- 
tic of Class I solutions towards the centre. Starting from 
equations 1221 and 1231 1. we can construct these new Class 
HI MHD solutions. In order to assure the existence of this 
Class HI solutions, A should be greater than 4 (Wang & Lou 
2006). As examples of illustration, we take K = 1.0, A = 5.0, 
H = [2- \+ [X^ - 4A)^/2]/2 and x, = 1.2, 1.5, 1.7, 2.2, 
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Table 4. Semi-complete Class I and II similarity twin MHD shock solutions matched 
with different outer magnetosonic critical point a;»(2). 



x.(2) 



Description 



^s{l) 



Xsii) 



0.90 



a;*(l) = 1.44 X 10"'' 0.8352 

mo = 2.8873 x lO-"* 

D = 3059 0.8297 



0.95 



1.00 



x,{l) = 1.88 X 10"* 0.7901 

mo = 3.7666 x lO""* 

D = 2344 0.7854 



V 



1.5552 -0.3852 



1.00 


1.6398 


-0.3034 


1.10 


1.8546 


-0.1103 


1.15 


1.9868 





1.20 


2.1401 


0.1198 


1.25 


2.3238 


0.2553 


1.05 


1.8414 


-0.1298 


1.10 


1.9615 


-0.0289 


1.15 


2.1070 


0.0869 



Columns 1 to 6 give the values of x* (2) where the Class I and Class II twin shock solutions 
cross the sonic critical line the second time, x*(l) and mo for Class I solutions, D for 
Class II solutions, the first shock location Xs{l), the second shock location Xs{2), the 
mass density parameter A and the speed parameter V. 
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Figure 21. Class I MHD shock solutions with A = 0.1 for 
a;»(2) = 0.9 and 1.0 for the ratio of the Alfven wave speed and the 
isothermal sound speed v^/a versus x for the solutions shown in 
Fig. 1201 The two corresponding values of a;»(l) are 1.44x 10""* and 
1.88 X 10~*, respectively. The corresponding inner MHD shock is 
located at Xs = 0.8352 and Xs = 0.7901, respectively. 



Figure 22. Class II MHD shock solutions with A = 0.1 for 
Xt = 0.9 and 1.0 in terms of —v{x) versus x by the heavy and light 
solid curves, respectively. The corresponding values of (D, Xs) are 
(3059, 0.8297) and (2344, 0.7854), respectively. The dash-dotted 
line represents the magnetosonic critical curve. The speed ratio 
v^/a versus x is shown in Fig. 1231 



respectively. In Figs. 1^ and 1291 we present these Class III 
MHD shock solutions. 



4.6 Two- Temperature Similarity MHD Shocks 

We have only studied so far the case of r = 1 in the preced- 
ing sections. This means that both the downstream and up- 
stream magnetofluids across an MHD shock have the same 
thermal temperature, i.e., the entire magnetofluid has same 
thermal sound speed Oi, = a^. In some astrophysical sys- 
tems (e.g., HIT regions around luminous massive OB stars), 
the factor of determining the ionization temperature, mean 
atomic mass, etc., would vary from an interior region to an 
exterior region. In this section, we shall discuss a magnetized 
gas medium of two different temperatures across an MHD 
shock. In our following consideration, the downstream sound 
speed Ud is thus no longer equal to the upstream sound speed 



ttu- In most astrophysical systems, the downstream sound 
speed should be faster than the upstream sound speed, while 
the downstream and upstream fluids can be regarded as 
isothermal fluids separately. With this in mind, we assume 
dd/ciu > 1 in our following MHD model analysis. 

Let us begin with the MHD shock solution for the case 
oi V = at X —> -l-oo. Here, we take Class I MHD solutions 
with a magnetosonic critical point a:,(l) = 0.103 and Class 
II MHD solutions with a central reduced density parameter 
D = 0.1 as examples of illustration. For different values of 
sound speed ratio Ud/au, we obtain the Class I and Class II 
MHD shock solutions through matching the upstream and 
downstream solutions in the speed-density phase diagram. 

In Fig. inland Fig. ED we display both classes of MHD 
shock solutions with their parameters summarized in Table 
5. In both classes of MHD shock solutions, as the ratio ad/a^ 
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Figure 23. Class II MHD shock solutions with A = 0.1 for 
Xt = 0.9 and 1.0 shown for the ratio of the Alfven wave speed 
to the isothermal sound speed v^/a versus x corresponding to 
the results of Fig. 1221 The corresponding values of (D, Xs) are 
(3059, 0.8297) and (2344, 0.7854), respectively. 
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Figure 24. Class I (solid curve) and Class II (dotted curve) 
twin MHD shock solutions with A = 0.1 for z» = 0.9 are shown in 
terms of —v(x) versus x. The dash-dotted line represents the mag- 
netosonic critical curve. The outer MHD shock locations Xs(2) are 
located at 0.95, 1.00, 1.10, 1.15, 1.20, 1.25 to match with var- 
ious asymptotic MHD solutions at a; — > -|-oo. The corresponding 
speed ratios vaJo- versus x are shown in Fig. 1251 



Figure 25. Class I (solid curve) and Class II (dotted curve) 
twin MHD shock solutions with A = 0.1 for a;* = 0.9 are 
shown in terms of the ratio of the Alfven wave speed to the 
isothermal sound speed va/o. versus x corresponding to re- 
sults of Fig. 1241 The outer MHD shocks locations is (2) are at 
0.95, 1.00, 1.10, 1.15, 1.20, 1.25 respectively to match with var- 
ious asymptotic MHD solutions at a; — > -|-oo. 
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Figure 26. Class I (solid curve) and Class II (dotted curve) 
twin MHD shock solutions with A = 0.1 and the magnetosonic 
point x, = 1.0 shown in terms of —v(x) versus x. The dash- 
dotted line represents the magnetosonic critical curve. The outer 
MHD shock locations a;s(2) are at 1.05, 1.10, 1.15, respectively 
to match with various asymptotic MHD solutions at x — > -|-oo. 
The corresponding speed ratio va/o, versus x is shown in Fig. 1271 



increases, the downstream shock location x^d, namely, the 
Mach number of the downstream shock becomes smaller. 
The mass parameter A increases with an increasing ad/a^ 
ratio. When the mass parameter A is less than 2, the up- 
stream solution is an MHD outflow breeze. The MHD breeze 
is stronger with a smaller value of A. In the case of the mass 
parameter A larger than 2, the inflow speed becomes weaker 
when A becomes smaller. 

We now investigate the case for the velocity parameter 
V ^ Q. For instance, we choose ad/au = 1.5 to study Class 
II MHD shock solutions. Different locations of x^d lead to 
different MHD shock solutions. In Fig. El and Fig. ESI we 



present different MHD shock solutions with the reduced cen- 
tral density parameter D = 10~*. 



5 NOTES AND DISCUSSION 

In this paper, we investigated quasi-spherical self-similar 
shock solutions for magnetized self-gravitating isothermal 
fluids in the semi-complete solution space. The relevant sim- 
ilarity MHD shock solutions are obtained and classified by 
comparing with resii l ts of previous analyses (iTsa i fc Hsul 
ll995l:IShu et alJl2002l: ILou fc ShenllJOoi IShen fc Loull2004l : 
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Table 5. Two-temperature Class I and II similarity MHD shock breeze solutions with 
different values of sound speed ratio ad/ciu- 



Parameter 


a-dlo-u 


A 


^sd 


^STi 


Vd 


Q-d 


Vu 


Oil 


a;.(l) = 0.103 


1 


1.8530 


1.5018 


1.5018 


0.5805 


1.2907 


0.1507 


0.8802 




1.3 


2.7594 


1.4834 


1.9284 


0.5750 


1.3150 


-0.4102 


0.6640 




1.5 


3.5226 


1.4612 


2.1918 


0.5682 


1.3452 


-0.6657 


0.6305 




1.8 


4.9048 


1.4414 


2.5945 


0.5622 


1.3729 


-0.9916 


0.6059 


D = 0.1 


1 


0.2976 


2.3960 


2.3960 


1.7176 


0.1563 


0.8273 


0.0676 




1.3 


0.4226 


2.3120 


3.0056 


1.6517 


0.1521 


0.5082 


0.0523 




1.5 


0.5288 


2.2786 


3.4179 


1.6256 


0.1505 


0.3913 


0.0487 




1.8 


0.7206 


2.2480 


4.0464 


1.6018 


0.1491 


0.2675 


0.0459 



Columns 1 to 9 contains relevant parameters of MHD shock solutions: a;»(l) for the 
magnetosonic critical point of Class I MHD solutions and D for the central reduced 
density parameter of Class II MHD solutions; the sound speed ratio of the downstream 
ad to the upstream a^ ; the upstream mass density parameter A\ x^d is the shock location 
in terms of downstream a^j; Xsu is the shock location in terms of upstream a^j; the 
downstream reduced speed Vd\ the downstream reduced density ad\ the upstream reduced 
speed Vu\ the upstream reduced density a^. 

Table 6. Two-temperature Class II MHD shock solutions for sound speed ratio a^/a^ = 
1.5 with D = 10"". 
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■"d 
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Vu 


Olu 


A 
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2.0 


3.00 


1.4407 


1.5753 X 10-6 


-0.2092 


4.1184 X lO-'^ 


3.3637 X 10-6 


-0.7544 


2.5 


3.75 


1.8608 


2.0284 X 10-" 


0.7814 


6.5518 X lO-'^ 


8.3927 X 10-6 


0.2820 


3.0 


4.50 


2.3032 


2.7282 X 10-'' 


1.6533 


1.0016 X lO-'' 


1.8592 X 10-5 


1.1984 



Columns 1 to 8 summarize relevant parameters for MHD shock solutions, x^d is the 
shock location in terms of downstream a^j; Xsu is the shock location in terms of upstream 
au', the downstream reduced velocity is Vd; the downstream reduced density is a^; the 
upstream reduced velocity is «„; the upstream reduced density is au; A is for the mass 
density parameter of upstream and V is the upstream speed parameter. 
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X (1)=0.7854 
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Figure 27. Class I (solid curve) and Class II (dotted curve) twin 
MHD shock solutions with A = 0.1 and a magnetosonic critical 
point x* = 1.0 shown in terms of the speed ratio Vj^/a versus 
X corresponding to the MHD shock results of Fig. 1261 The outer 
MHD shock locations Xg{2) are at 1.05, 1.10, 1.15, respectively, 
to match with various asymptotic MHD solutions at x — > +oo. 



iBian fc Lotl 1200,4 IYu fc Loull200,4 l. For a closely relevant 
similarity analysis on a magnetized polytropic gas, the 
reader is referred to Wang & Lou (2006). Such mEECC sim- 



ilarity solutions, which depict an expanding envelope with a 
concurrent collapsing core, are derived from nonlinear MHD 
equations and may be applied to various astrophysical MHD 
processes with appropriate adaptations. 

As noted by Yu & Lou (2005), the velocity and den- 
sity profiles of an EECC shock model agree better with 
observations of cloud B335 than the results of the EWCS 
model. The magnetized EECC shock solutions qualitatively 
retain this essential feature of nonmagnetized EECC solu- 
tions (Shen & Lou 2004). The inclusion of a random mag- 
netic field is more realistic and can be further utilized to 
model diagnostic features such as synchrotron emissions etc. 
Due to current observational limits, no empirical magnetic 
field information is immediately available for the most in- 
ner core of cloud system B335 (e.g., Wolf et al. 2003). Our 
MHD EECC shock solutions may give a rough estimate for 
the magnetic field strength in the most inner core of star 
forming cloud B335 as an example. For the starless cloud 
system B335, we may take an infalling region of a spatial 
scale ~ 1.5 x 10'* AU, a cloud mass of ~ 4:Mq, a random 
magnetic field strength of ~ 134/iG (e.g.. Wolf et al. 2003) 
and a thermal sound speed a ~ 0.23 km s"'^. For a quasi- 
spherical accretion MHD shock to exist within such a cloud 
system, the core magnetic field strength may reach ~ ImG, 
which is several times strong than our previous estimate 
of ~ 300/iG, depending on the strength of such an MHD 
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Figure 28. Semi-complete Class III MHD sliocic solutions in 
terms of —v{x) versus x with A = 5.0 and K = 1.0 for shock 
location Xs = 1.2, 1.5, 1.7, 2.2, respectively. Asymptotic MHD 
solutions (iv) given by equations 1221 and I2ij| are specified near 
the origin. The dash-dotted line represents the magnetosonic crit- 
ical curve. The speed ratio va/h versus x is shown in Fig. 1291 



Figure 30. Class I (light solid curves) two-temperature MHD 
shock breeze solutions with A = 0.1, x* = 0.103 and differ- 
ent density ratio a^/ciu ^ 1- Class II (heavy solid curves) two- 
temperature MHD shock breeze solutions with A = 0.1, D = 0.1 
and different density ratio a^/a^ ^ 1. The dash-dotted line rep- 
resents the magnetosonic critical curve. The reduced transverse 
magnetic field strength b{x) versus x is shown in Fig. 1311 




^0.4 




Figure 29. The ratio of the Alfvcn speed to the isothermal sound 
speed va/ci versus x with A = 5.0 and K = 1.0, corresponding 
to the MHD shock results in Fig. I^U Class III MHD shock so- 
lutions with the MHD shock location at Xs = 1.2, 1.5, 1.7, 2.2, 
respectively. 



shock. Observations indicate that some T Tauri stars may 
have surface magnetic field strengths of the order of lO^G, 
which is much greater than our estimate. Possibly, our quasi- 
spherical solutions may not be valid when the spatial scale 
involved becomes much smaller than ~ 100 AU; within such 
a scale, a circumstellar disc as well as the stellar dynamo 
process would dominate among various processes to amplify 
the protostellar magnetic fields (n.b., a small-scale circum- 
stellar disc and a relevant bipolar outflow would be regarded 
as additional features in our quasi-spherical mEECC shock 
scenario) . 

Such mEECC similarity solutions with a random mag- 
netic field may also be applied to the asymptotic giant 
branch (AGB) phase or post-AGB phase in the late evo- 
lution stage of a low-mass main-sequence star before the 



Figure 31. Reduced transverse magnetic field strength b(x) 
versus x corresponding to the results in Fig. 1301 Class I two- 
temperature MHD shock breeze solutions with A = 0.1, Xt = 
0.103 and different density ratio a^/oLu s? 1. Class II two- 
temperature MHD shock breeze solutions with A = 0.1, D = 0.1 
and different density ratio a^/au JS 1. 



gradual emergence of a planetary nebula (PN) system with 
a central magnetized white dwarf of a high surface temper- 
ature in the range of 1 x 10^ ~ 2 x lO^K. The timescale of 
this evolution 'gap' between these two phases is estimated 
to be ~ lO'^yrs. Planetary nebulae and proto-planetary neb- 
ulae (pPNs) are believed to be the ultimate evolutionary 
stages of low- and intermediate-mass stars (M^ SMq). PNs 
and pPNs appear on the sky as expanding plasma clouds 
surrounding a luminous central star. Such expanding clouds 
have a typical asymptotic speed of ~ 10 — 20 km s^^; the 
relevant mass loss rates fall within the range of ~ 10~** to 
~ 10~ Mq yr"'^. With an insufficient nuclear fuel supply 
from a certain epoch on, the central region starts to contract 
and collapse while the outer envelope continues to expand 
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Figure 32. Class II MHD shock solutions with A = 0.1 and 
D = 10~^ shown in terms of —v{x) versus x by light solid curves. 
The dash-dotted line represents the magnetosonic critical curve. 
The reduced transverse magnetic field b(x) versus x is shown in 

Fig.oa 



_Ql-5 




Figure 33. The reduced transverse magnetic field h{x) versus x 
corresponding to the Class II MHD shock solutions with A = 0.1 
and D = W^^ shown in Fig. El 



into a massive wind; it is highly likely that the central col- 
lapse is accompanied by an outward energetic surge to chase 
the slowly expanding envelope. Thus a slow, dense stellar 
wind expelled during the AGB phase is followed by a fast, 
tenuous magnetized wind driven off the collapsing proto- 
white dwarf during the PN phase. The collision of magne- 
tized winds with two different flow speeds would inevitably 
generate an MHD shock. In this process, the outer expansion 
removes stellar envelope mass companied by a proto white 
dwarf produced at the centre by central infall and collapse. 
Sufficiently far away from the initial and boundary condi- 
tions, the system may gradually evolve into an MHD phase 
representable by an MHD EECC similarity solution during 
a timescale of a few hundred to several thousand years. As 
noted out by Bian & Lou (2005), the dimensionless mass 
accretion rate mo should be in the range of ~ 10""^ — 10~*. 
In a PN, the dimensionless magnetic parameter A is about 
0.003. Our Class I twin MHD shock solutions appears to 



have the desired reduced mass accretion rate towards the 
collapsed core. 

Observations grossly indicate a magnetic field strength 
of ~ IG near the stellar surface at r ~ lAU during the AGB 
phase, decreasing approximately to ~ 10~^G at r ~ 10"^ AU. 
In the quasi-spherical self-similar MHD expansion regime of 
our model analysis, the transverse magnetic field strength 
scales as B±^ oc r~^ (e.g., Lou 1993, 1994); for an expanding 
MHD shock travelling outwards, the magnetic field down- 
stream would be enhanced by the existence of such an MHD 
shock. Our model estimates are roughly consistent with 
the measured magnetic field variations observed for AGB 
stars. The magnetic field strength associated with an MHD 
EECC solution in the innermost collapse region scales as 



B 



-1/2 



If we take the strength of a surface magnetic 



field as ~ IG at r ~ lAU and the radius of a proto-white 
dwarf as ~ 6000 km, then the magnetic field strength es- 
timated at the surface of a proto-white dwarf is roughly 
lO'^G. An accretion MHD shock would further increase this 
inner downstream magnetic field. This enhancement is de- 
termined by the MHD shock strength. Much stronger surface 
magnetic fields of a magnetic white dwarf (~ 10^ — 10^ G) 
might be generated and sustained by intrinsic stellar MHD 
dynamo processes driven by the convective differential rota- 
tion inside an AGB star (e.g., Blackman et al. 2001). 

We therefore hypothesize that dynamical evolution of 
an mEECC phase of around or less than a few thousand 
years may be the missing linkage between the AGB or post- 
AGB phase and the gradual appearance of a PN. Depending 
on physical parameters of the low-mass progenitor star, it 
may also happen that the degenerate CO core collapse and 
subsequent material infalls during an mEECC phase lead to 
an eventual core mass exceeding the Chandrasekhar mass 
limit of 1.4Mq and thus induce a supernova with an intensity 
determined by the actual central mass accretion rate. 

For the Crab Nebula as an example of magnetized SNR, 
we model the nebular magnetic field as Bx oc r~^ during the 
outer envelope expansion phase, ignoring complex interac- 
tions between the central magnetized relativistic pulsar wind 
and the inner nebula. In the presence of shocks, as a fiow 
passes through a shock from downstream to upstream, the 
magnetic field strength would decrease by a factor of a few. 
Thus, the magnetic field strength in the fiow would change in 
the presence of shocks as compared with mEECC solutions 
without shocks, although not significantly. For a neutron 
star of a radius ~ 10 km and a typical dipolar magnetic field 
of strength ~ 10^° G (we exclude for the moment those un- 
usual high-field magnetars, such as SGRs and AXPs whose 
magnetic field strength may reach as high as 10^* ~ 10^^ G), 
the magnetic field would decrease to ~ 10~"^G at a distance 
of about several parsecs (Yu & Lou 2005). When an MHD 
shock is included and for the same neutron star magnetic 
field, the outer nebular magnetic field would become weaker 
~ 10~^G at a distance of about several parsecs. These es- 
timates more or less agree with the magnetized envelope 
expansion portion of our mEECC shock solutions. In addi- 
tion to this order of magnitude agreement with observations 
for magnetic field strengths, mysterious central structures 
of the Gran Nebula like wisps and knots are likely produced 
by reverse fast MHD shock waves as a result of slightly in- 
homogeneous pulsar wind streams emanating from the fast 
spinning pulsar (e.g., Lou 1998). These quasi-stationary re- 
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verse MHD shocks in space can effectively produce relativis- 
tic particles and synchrotron emissions from relativistic elec- 
trons gyrating rapidly around magnetic field lines surround- 
ing these compact objects. The familiar triple ring structure 
of SN1987A may be associated with MHD shocks (Tanaka 
& Washimi 2002). 

Unlike the Crab Nebula, which is primarily powered by 
the rotational energy through a relativistic pulsar wind, the 
slowly rotating magnetars are thought to be powered by the 
decay of extremely intense magnetic field. The light curves 
may thus represent an adiabatically expanding population 
of electrons accelerated at a particularly active phase, which 
could be produced by the ejecta colliding with a pre-existing 
shell. Such a shell is naturally made by SGR itself, because 
its quiescent wind of a luminosity --^ 10"''' erg s~^ will sweep 
up a bow shock of a stand-off distance ~ 10^^ cm as it moves 
through the interstellar medium (ISM) at a typical neutron 
star velocity of ~ 200 km s~^. If this pre-existing shell is hit 
by giant flares from a SGR with an energy of ~ lO**^ — 10'*'' 
erg, MHD shocks will naturally occur and be swept outward, 
resulting in a violent episode of relativistic particle acceler- 
ation that deposits much of the energy into a steadily ex- 
panding synchrotron-emitting shell after giant flares. Polar- 
ized emissions from SGRs also indicate that magnetic field 
should play a crucial role in understanding SGR phenomena. 

As massive OB stars turn on their core nuclear reactions 
in magnetized molecular clouds, the interstellar medium of 
surrounding neutral hydrogen gas strongly absorbs ultravio- 
let radiations from OB stars and becomes ionized to form lu- 
minous HII regions (e.g., Stromgren 1939; Osterbrock 1989; 
Shu et al. 2002). As the ionization front sweeps through the 
neutral cloud, magnetized gas flows driven by pressure gra- 
dients develop between HII and HI regions as well as within 
HII regions, and MHD shocks can naturally emerge in these 
processes. 

For astrophsical systems of much larger scales such as 
clusters of galaxies (e.g., Sarazin 1988; Fabian 1994), the 
similarity MHD shock solutions may be valuable for under- 
standing a certain phase of their evolution involving mag- 
netic fields (e.g., Hu & Lou 2004; Lou 2005). Chandra obser- 
vations show that the continuous blowing of bubbles by the 
central radio source would lead to the propagation of shocks 
seen as the observed fronts and ripples, gives a rate of work- 
ing which balances the radiative cooling within the cluster 
core (Fabian et al. 2003). With appropriate adaptation, our 
MHD shock model may offer interesting interpretations of 
such phenomena. 
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APPENDIX A: DERIVATION OF MHD 
PERPENDICULAR SHOCK CONDITION 

We have the following basic definitions and relations 

pi ui X pi Bi a 

Momentum conservation 1301 can be written as 
Xpi + X^Bf/{8n) + piuj/X =pi + B!/{8n) + pmi . 
Using the following two relations 

8n /3i ' 

the foregoing momentum equation becomes 
Xpi + X^pi/Pi + piMl/X = pi + pi//3i + piMl , 
or simply 

X + A:V/3i + M^/X = 1 -H l//3i + Ml . 
It then follows immediately that 
(X-l) + {X- 1){X + l)//3i = Ml{X ~ l)/X . 



2 
2 2 ^1 71 *-2 

piUi = pia — r- = piMi , 
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The degenerate case of X— 1 = would be trivial. Removing 
this factor X — 1, we readily arrive at 

1 + (X + i)//3i = m!/x . 

This is the quadratic equation governing the density ratio 
X, exactly the same as equation i'SZl in the main text. 
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